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Abstract 

We use central differences to solve the time depen- 
dent Euler equations. The schemes are all advanced 
using a Runge-Kutta formula in time. Near shocks 
a second difference is added as an artificial viscos- 
ity. This reduces the scheme to a first order upwind 
scheme at shocks. The switch that is used guaran- 
tees that the scheme is locally TVD. For steady state 
problems it is usually advantageous to relax this con- 
dition. Then small oscillations do not activate the 
switches and the convergence to a steady state is im- 
proved. To sharpen the shocks different coefficients 
are needed for different equations and so a matrix 
valued dissipation is introduced and compared with 
the scalar viscosity. The connection between this ar- 
tificial viscosity and flux limiters is shown. Any flux 
limiter can be used as the basis of a shock detector for 
an artificial viscosity. We compare the use of the van 
Leer, van Albada. minmod. superbee and the ''aver- 
age" flux limiters for this central difference scheme. 
For time dependent problems we need to use a small 
enough time step so that the CFL was less than one 
even though the scheme was linearly stable for larger 
time steps. Using a TVB Runge-Kutta scheme yields 
minor improvements in the accuracy. 
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1 Basic Scheme 

The basic elements of the scalar dissipation model 
considered in this paper were first introduced by 
Jameson. Schmidt, and Turkel [2] using an explicit 
Runge-Kutta time integration scheme. The space dis- 
cretization is based on central differences with an ad- 
ditional artificial viscosity. In this section the basic 
scheme is briefly reviewed. 

Consider the Euler equations in the form 

W t + f 2 = 0. (1) 

where W is the three-component vector of conserved 
variables, and / is the flux vector. The independent 
variables are time t and Cartesian coordinate i. In a 
cell-centered, finite- volume method, (1) is integrated 
over an elemental volume in the discretized compu- 
tational domain. Equation (1) can also be written 
as 

W t + AW X = 0 

where A is the flux Jacobian matrix defined by A - 
df/dW . 

To advance the scheme in time we use a multistage 
scheme. A typical step of a Runge-Kutta approxima- 
tion to (1) is 

W^ = W {0) -a k ^ -AV], (2) 

where D is the spatial differencing operator, and 
AV represents the artificial dissipation terms. The 
derivatives of the fluxes are approximated by central 
differences. In the form presented here the scheme 
can not have greater than second order accuracy in 
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r.iino for nonlinear problems. For steady state prob- 
lems the time accuracy is irrelevant and the form of 
(2) requires only two levels of storage. If one wishes 
to obtain higher accuracy in time for nonlinear prob- 
lems. then one can use any formula from standard 
numerical ODE theory. In particular the classical 
Runge-Kutta scheme will give fourth order accuracy 
using four stages but will require more storage than 
(2). We will also consider Runge-Kutta forms that 
preserve the TVD nature of the scheme when the spa- 
tial operator is TVD [4]. [5]. In all cases, the spatial 
accuracy is determined only by the accuracy of the 
operator D to the derivative. 

The dissipation terms are a blending of second and 
fourth differences. That is. 

AV = (D 2 - D 4 ) W, 

where 

u 2 w = v[(A, +i € ( .;^)A]^. (3) 

= v[(a, + , € '1 ) ,) AVA] Wi. (4) 

and A . V are the standard forward and backward 
difference operators respectively. The variable scaling 
factor A is chosen as 


2 Shock detectors and flux 
limiters 

In order to sec the effect of v we first define 

4>i = 1 - v x . (7) 

As shown in [6] </> can be interperted as a flux limiter, 
though its properties for central difference schemes is 
slightly different than for upwind schemes. The value 
of <f> is usually taken as a function of r where 

A_ , , 

T = = -t— I®) 

Ui+i - Ui A+ 

According to the TVD theory for a scalar equation in 
one dimension the artificial viscosity can sometimes 
be negative , see (20). However, for multidimensional 
vector equations with central differences we prefer 
to be conservative and choose the artificial viscosity, 
e* 2 ), to be positive and so we set 

^ = ( 9 ) 

For the fluid dynamic equations we choose the pres- 
sure as a representative of u. The artificial viscosity 
used in the original algorithm was 


= 9 l^i + ^i+l] * (5) 

where A* is proportional to the wave speed. The co- 
efficients e {2) and e u> are adapted to the flow and are 
defined as follows: 

= n- 2) max(t/j.i/i+i). (C) 


e'+i = max [°‘ (* !41 - • 

where k 12} and k ( 4} are constants to be specified. 

The parameter v is a shock detector. We shall an- 
alyze ways of defining v in detail in the next section. 
The purpose of this second difference viscosity is to 
introduce an entropy-like condition and to suppress 
oscillations in the neighborhood of shocks. Ideally the 
value of v should be one at shocks and be negligible 
in smooth regions of the flow. The fourth-difference 
dissipation term is basically linear and is included to 
damp high-frequency modes and allow the scheme to 
approach a steady state. Only this term affects the 
linear stability of the scheme. Near shocks it is re- 
duced to zero. For time dependent, flows, the fourth 
order dissipation is not very important and k {4) will 
usually be small or zero. 


Pt+l 2 Pi 4 pi — 1 

1 4- 2 Pi 4* pt-i 


( 10 ) 


and zq+i /2 = max{vi,Vi+x). We note that with this 
definition of v that <j> is not a function of r. We shall 
demonstrate in the result section that this switch 
gives rise to oscillations in the flow field. 

In order to connect this artificial viscosity with flux 
limiters we first consider the van Leer flux limiter 
given by 




r 4 |r| + £ 
14 |r| -F € 


( 11 ) 


Where e is added to prevent the switch from being 
activated by noise. This € is mainly needed for steady 
state calculations. Then after multiplying (11) by 
|A+| we get 


1 -<f>r(T) = 


(A + — A-)sgn(A + ) 
|A_| + |A + | + e|A + r 


( 12 ) 


Reverting back to the notation of pressure and mod 
ifying the € term we get 


^ = |i - <M r )l = 


lgr+1 - 2 p, +p,-i| 

|p.+i — pU + b. -p.-i| + e' 


For dimensional consistency we wish to choose e to 
depend on the pressure. So we choose e — e(pi+i + 

~pi + p.-i) • 
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Hence. 

v _ jP»+i - 2p, +p,-i| 

|p«+i - Pi | + I/'. - 7^i — 1 1 + «(p<4- 1 + 2j)i +p,-i) 

(13) 

There is no special need to base the artificial viscos- 
ity on the van Leer flux limiter. It is just coincidental 
that the resultant viscosity v closely resembles the 
original artificial viscosity. (10). Another alternative 
is the van Albada flux limiter. 


We note that this limiter approaches zero for large 
values of r. while most limiters approach 2 as r in- 
creases. Using a similar derivation we find that the 
artificial viscosity associated with the van Albada 
limiter is given by 


any difference in (8) if r is a forward difference over 
a backward difference or a backward difference over 
a forward difference. For smoothness we now want 
0 (1) = 0. Of the above limiters only the first version 
of Van Albada and superbee have this property. It 
follows from the analysis of [1] that an upwind scheme 
can be considered as a symmetric interpolation fol- 
lowed by a upwind convection operator. A central dif- 
ference scheme can be represented as a downwind in- 
terpolation followed by a compensating upwind con- 
vection and so the total operation is symmetric. 

3 The TVD Property 

Consider the one-dimensional scalar conservation law 
|[u(i,*)] + ^[/(«(*,*))] = 0, (17) 


(p«+ 1 ~ 2p, +p,_i) 2 

(Pi+i - Pi ) 2 + (p, -p,-i) 2 + e(p<+ 1 + 2 pi+pi-\ 

(15) 

There is a second version of the van Albada flux 
limiter used in the literature. 

0 < r < 1. (16) 

1 H- r 

Other limiters used are minraod 

<t>i[r) — max(min{r , 1). 0), 
and superbee 

<pi(r) = maz(min(2r. 1). mm(r, 2). 0), 

see [7]. 

We shall also consider the "average" flux limiter 

<f>i(r) = minmod({l + r)/2,2mmmo(!(l.r)). 

For each of these limiters there is a corresponding 
•artificial* viscosity. 

For an upwind flux limiter we have 0(£) = ;0(r) . 
Kuynh [1] has shown that the resultant scheme is 
second order if 0(1) = 1/2. However, for a central 
difference scheme we have 

0 central — 1 “ t 7 ; — 1 [1 01 1 

and so 

_ f <pi if 0 < 1 . r < 1 

tpcentral - <y , _ Q. if 0 > 1 . r > 1 


where 

i -oc < x < oc, t > 0. 

Let v{t) = {vj(t)} be the approximate solution of (17) 
and consider the semidiscrete equation 18 

< 18 > 

with 

A i> i+ i = (Av) i+ i = Vi+i(«) - Vi(t). 

A 3 is a third-difference operator defined as 

A 3 v i+ ^ = v i+2 {t) - 3v i+ i(t) + 3vi{t) - Vi-i(t). 

The terms on the right hand side of (18) represent 
second- and fourth-difference numerical dissipation 
terms, with k ^ a constant. Define 

s i+ \ =sg n ( Al, <+0- 

where sgn represents the signum function. We first 
shift the indices by one in (18) and subtract (18) from 
the resulting equation. We then multiply the result 
by and sum over all i. Noting that = ±1, 

so 5*^ = 1. and 

S, + 1 Ar, + i = | Av i+ i . 

Let TV denotes the total variation as given by 


With the van Leer and first Van Albada flux lim- TV = V* At» <+ J . 

iters one finds that 0(i) = <f>(r). i.e. it. doesn't make < 


3 



inii mum u mi 


we then obtain 


i£i Al ’^i = 


2Ax 2. S i+i ( S i~i S »+i) Av, + j | A "*+i 
■^2Ax^t 5 «+i ( S > + 3 — ^ s t+4 + S >-$) AV *+a | 

-CE,( s H‘ 2 M +s h) i? i+ iA 3 t; i+ ^ 

We stress that the last term will not help for TVD. 
Its purpose is to eliminate high frequencies and ac- 
celerate convergence to a steady state. Hence, we 
want this contribution to be zero. This can be ac- 
complished if we demand either 

{ s l+ 3 -2s, + i +s t _ i =0 

°l+l = 0, 

We are then left with 

it {TV) = (19) 

a/, + i I 

_ 2^i -S i+^( S *+? _S >-pAv i+ ^ | A,; ‘+i 

-*i+$ ( s i+3 - 2s,+^ + )Q i+ i |Ai\ + ^ | 

Thus, a sufficient condition that the total variation 
not increase is that each term in the summation of 
(19) must be positive. This means that the scheme 
is TVD if 

- s, + i (s.+jj - 2s i+ i + s,_i) Q i+ i > (20) 

( \ 

This is the inequality obtained in [6]. 

When driving the solution to a steady state one 
frequently finds that it is not advantageous for the 
scheme to be TVD. The reason is. that with TVD 
schemes the switches are frequently being turned on 
and off due to local noise. For steady state calcula- 
tions this causes the convergence to halt at some error 
level and a limit cycle results in which the residual 
oscillates about some level instead of decreasing. To 
prevent this from occurring we wish to prevent the 
switch from being activated for small oscillations or 
small discontinuities. The inequality (20) was ob- 
tained by demanding that the solution be TVD and 
so each term on the right hand side of (10) was neg- 

ative independent of the size of - ■ . Instead we 
shall only demand that the solution be total variation 


bounded (TVB). Now, each term on the right hand 
side of (19) can be positive as long as it is bounded 

by a constant times At/ l+ i| . Since s^+i /2 Is equal 
to plus or minus one we want 


A/i+i 


- Qi+1/2 + 

a a positive constant. 

We shall choose 


Qi+1/2 = v i+ 1/2 , - ( 22 ) 

£ * V i+\ 

This is similar to (5,6) with k (2) = 1/2, */,-+ 1/2 = 
max(vi, Vi+i), and A = ^ ^ • We then rewrite (22) 


(1 - "i+ 1 / 2 )- 


This implies that if is small then we don’t 

need to turn on the artificial viscosity parameter v. 
Only when — -+ 00 do we need that v -* 1 . In 


Only when 


A/ i 

a steady state the shock speed, is zero and so 

(23) is satisfied for any positive Q. So for a steady 
state scalar equation the TVD property is trivial. 

Hence, for a steady state problem we do not need the 
complete TVD theory. However, in this case the solu- 
tion is also trivial. Moreover, the theory for systems 
is still not adequate for our purposes. Alternatively, 
we choose v to depend on the strength of the shock. 

£\v. For weak shocks Av is small and we can choose 
£ near one. For strong shocks Av is large and we 
want £ to be small so that v is a TVD switch. For 
the fluid dynamic equations we replace the vector v 
by the scalar pressure, p. 

To find such a v we use (13) When e = 0 we get 
the TVD switch (13) while with e = 1 we obtain a 
perturbation of the original switch, (10) for transonic 
flows. This switch treats the two sides of the shock 
asymmetrically depending on whether p* is to the left 
or right of the shock. Thus, we replace it by 

_ \pi+l ~~ tyx -bPi-l[ ____ 

V% ” |p»+l ~Pi | + | Pi ~P*-l| + c 77iax(pi+i,p*,pi-i) 

In practice the switch that we use is 

jgt+i ~ 2pj + Pi-il 

1/1 ” (1 - e)(|pi+i - Pi | V | Pi — Pi-i I) + + 2 P< + P*-i) 

(24) 
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Wo wish to choose t automatically based on the 
shock strength. One possibility for e is 

\max(pi-2*P»-i*P».Pt+i-P»-2) 



where a is a free parameter. A reasonable range is 
(7 = 1/2 to a = 1. 

For small oscillations doesn't vary much and so 
€ is slightly less than 1 . For large oscillations e is 
equal to the relative jump across the shock. Consider 
a perfectly resolved discontinuity going from pi to pr 
with pr < Pl - Then. 


Pl ~ PR 

(Pl -Pr) + *Pl Pl 


Let. 


Q = 


PR 
Pl 

Combining these we find that. 


1 


1 + 


1-9 


Hence, for small q (i.e. for large discontinuities) v ~ 
1 — q a . Hence for both very weak shocks and very 
strong shocks the left hand side of ( 23 ) is small i.e. 
(1 - i/) At; -* 0 as Av -» 0 and also as Av oc . 
This discussion has concentrated on the theoretical 
basis of the algorithm. In practice the formula ( 24 ) 
is used for transonic, supersonic and hypersonic flow 
regimes. 

In this section we have written the flux limiters and 
artificial viscosity in terms of the pressure variable 
which is appropriate for inviscid fluid dynamics. In 
the next section we shall consider matrix viscosities. 
With a matrix viscosity one can base the limiter in 
each characteristic field on a different quantity. 


4 Matrix Viscosity 

In the above discussion we have discussed a scalar 
equation. In the original algorithm, this procedure 
was applied to each equation with the same e. The 
coefficient A in ( 5 ) was chosen as equal to the spectral 
radius |u| -f c while v was the same switch that de- 
pended on the pressure, for all the equations. For 
time dependent flows this presents several difficul- 
ties. first as seen in the result section there is exces- 
sive smearing since the same coefficient is used for all 
waves and is proportional to the fastest wave speed. 
Second, pressure is continuous across a contact dis- 
continuity and so a pressure based switch will not 
sense a contact. We therefore replace the scalar dis- 
sipation with a matrix dissipation, i.e A in ( 3 . 4 ) is now 


a matrix valued function. We first define a function 
of a matrix A. We assume that A can be diagonal- 
ized so that TAT " 1 is diagonal. We then define our 
function as 


/(A) = A,/ + th. - A 3 )[^£i + £>] 

+^^-[E 3 + ( 1 -l)E,} 

where 



/ -u 1 
Ez = ( “U 2 u 

\ -uH H 



0 

0 

0 


0 

0 

0 


E< = 


0 

a 2 

a? 

2 


0 

-u 


0 

1 

u 


Whenever the matrix A can be diagonalized, i.e. 
D = T" l AT is diagonal then a function of the ma- 
trix is defined by f(A) = Tf(D)T ~ l , and f{D) is the 
function f applied to each element of the diagonal of 
D. Let the coefficients Ai,A 2 and A3 be functions of 
the eigenvalues of A. If Ai = u + c, A2 = u - c, A3 = u 
then we recover the matrix A. When the A’s are 
the absolute value of the eigenvalues we obtain the 
absolute value of the matrix A. In general, Ai,A 2 
and A3 should not be exactly equal to the eigenval- 
ues of A since at sonic points or stagnation points 
an eigenvalue is zero and hence no artificial viscosity 
would be added. Hence, the A’s have a lower limit 
of 0 . 2 |tt + c|. This procedure also allows one to se- 
lect different switches for each eigenvalue. In partic- 
ular we shall base the switch for the nonlinear fields, 
with speeds X x and A 2 on the pressure. However, the 
pressure is continuous across a contact discontinuity. 
Hence, the switch for the linear field. A3 is based on 
the temperature, T = * . Putting these options to- 
gether we choose the A’s equal to e (2) and e (4) times 
the limited absolute value of the eigenvalues, see( 3 , 4 ). 



5 Results 

The results were all obtained using a multistage 
Rimge-Kutta scheme (2) to advance the solution in 
time. For most of the computational results the origi- 
nal Runge-Kutta coefficients [2] were used, an = 1/4, 
ao = 1/3, az = 1/2. a 4 = 1. Shu [4] introduced an- 
other set of coefficients to guarantee that the scheme 
is TVD in time but is only first order accurate. The 
three stage scheme has coefficients, a\ = 1/9. a 2 = 
1/3. aa = 1 while the four stage scheme has coeffi- 
cients, Qi = 1/16, ao — 1/6. ctz — 3/8. a 4 = 1 . 
The more stages that are used the larger the time 
step allowed by stability requireme nts. However, we 
found that using larger time steps introduced oscilla- 
tions into the solution. In practice we chose CFL = 
.75, and so there was no advantage to using the four 
stage scheme. Shu [5] also introduced higher order 
schemes for time dependent equations that are still 
TVB. These schemes can no longer be written in the 
simple form of (2). Instead each stage requires the 
use of the dependent variables and fluxes at previ- 
ous stages and and so more information needs to be 
stored. 

We solve the one dimensional Euler equations in 
the domain 0 < x < 10. The initial conditions 
are u = Ora/s.T — 300K everywhere. The initial 
pressure is discontinous with a ratio of p = 20 for 
0 < x < 5 to 1 for 5 < x < 10. The density and 
total energy are then calculated from the ideal gas 
law with 7 = 1.4. 

We first consider the standard central difference al- 
gorithm with a scalar viscosity and the original switch 
(10) and the original Runge-Kutta coefficients with 
CFL = 0.75. The first figure is a plot of density as a 
function of x at a nondimensional time of 0.2 . Large 
oscillations appear both between the rarefaction wave 
and the contact discontinuity and between the con- 
tact and the shock. In figure 2 the density is plotted 
with the standard switch replaced by the van Leer 
based switch (9.11). The change in the switch has 
eliminated all oscillations since the scheme is TVD 
for the scalar case with this switch [6]. There are still 
some small oscillations in the rarefaction and the con- 
tact is very smeared. In figure 3 we show the same 
case using the matrix viscosity. The switch for the 
nonlinear waves is based on the pressure as before. 
Since the pressure is continuous across a contact dis- 
continuity the switch for the entropy wave is based 
on the temperature, though one could also use en- 
tropy. The smearing near the contact is considerably 
reduced but there remains a small glitch near the 
sonic point. Using the van Albada(l) based switch 
improves the treatment of the sonic point. The use 


of superbec for the nonlinear wave introduced new 
oscillations as seen in figure 4. We conclude that for 
the central difference schemes superbee should never 
be used for the nonlinear waves. The results with 
minmod was similar to the van Leer viscosity but 
with a slightly less sharp shock. In all cases the head 
of the rarefaction wave was smeared out. In figure 

5 we present the density when superbee is used for 
the linear wave while van Aibada(l) is used for the 
nonlinear waves. We also used these schemes with 
the e as given in (11). For the van Leer limiter we 
could choose e = .1 without significantly derading the 
results while for van Albada(l) we had to choose e 
about 0.005. For the steady problems we can use the 
van Leer limiter with e — .1 and still get monotone 
profiles. 

The cases presented until now were with the orig- 
inal four stage Runge-Kutta coefficients and CFL = 
0.75. Raising the CFL number introduced oscilla- 
tions. We next tried the first order scheme suggested 
by Shu [5] but still got oscillations when the CFL 
was larger than 1. We then used the third order 
Runge-Kutta scheme suggested by Shu [4]. Using 
the same switches for both the linear and nonlinear 
switches and these third order Runge-Kutta coeffi- 
cients resulted in a sharper profile but some oscilla- 
tions. Hence, in figure 6 we present the results for 
Shu*s third order scheme in time, using superbee for 
the linear field and the van Leer viscosity for the non- 
linear field. Figure 7 presents the same case as figure 

6 but with the CFL raised to 0.95. This introduced 
a small oscillation near the sonic point but otherwise 
was very satisfactory. It is interesting to note that 
with the scheme of Lerat and Sides [3] the solution 
becomes less oscillatory as the time step is increased. 
In our last case we consider the effect of using differ- 
ent variables for the switches. Until now the switch 
for the nonlinear fields has been based on the pressure 
while the switch for the linear field has been based on 
the temperature. We now plot the results when each 
characteristic field has an artificial viscosity switch 
based on that characteristic variable. In figure 8 the 
density is plotted for this case using Shu’s third or- 
der Runge-Kutta coefficients, CFL = 0.75, the super- 
bee limiter for the contact discontinuity based on the 
linearized entropy variable A p - c?Ap, and the van 
Leer limiter for the acoustic variables A p + pc An and 
A p - pcAu. We see that there is an additional im- 
provement in the resolution of the contact discontinu- 
ity. The solutions presented are all for the time t=0.2. 
Watching the time evolution one sees that there are 
many oscillations that occur in the initial breakup of 
the discontinuity into a shock and a contact. These 
oscillations disappear as the solution progresses. 
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We finally consider a steady state two dimensional 
calculation. We solve for the turbulent flow about a 
blunt cone using a Baldwiu-Lomax turbulence model 
at A/ x = 25.. a = 0 . The grid is 400 x 80. The 
geometry is shown in figure 9a. In figure 9b we plot 
the pressure and Mach number along the coordinate 
line directly in front of the cone. We choose e = 
0.9 in (24). With this high value of e there are only 
three points in the shock and no overshoots even at 
this hypersonic speed. If e is chosen equal to one 
(i.e. original switch (10)) the code does not converge. 
With smaller values of e the convergence is slowed. In 
figure 9c we show the convergence rate for this case 
on the three grids used in the FMG version of the 
multigrid. 

6 Conclusion 

• Using a central difference scheme with an artifi- 
cial viscosity we can duplicate most of the prop- 
erties of upwind TVD schemes. Solving the one 
dimensional time dependent Euler equations we 
get high resolution solutions for the shocks and 
the contact discontinuity. The main ingredients 
are an improved shock locator and a matrix arti- 
ficial viscosity. This shock locator can be based 
on any of the upwind flux limiters. Superbee is 
the best for the contact while either van Leer or 
van Albada(l) is best for the nonlinear waves. 
Further minor improvements can be obtained by 
using a high order TVD Runge-Kutta scheme 
and basing the switch on the characteristic vari- 
ables. The TVB Runge-Kutta schemes is slightly 
more accurate than the standard Runge-Kutta 
schemes. 

• There are major differences between the time de- 
pendent problem and the steady state problem. 
For the time dependent problem it was necessary 
for the scheme to be TVD-like in order to avoid 
oscillations. However, for the steady state prob- 
lem we use a coefficient for the artificial viscos- 
ity that is considerably below that required for 
the solution to be TVD and still get monotone 
shocks even for very strong discontinuities with 
about four points in the shock layer. Further- 
more. high order TVD schemes frequently slow 
down the convergence to the steady state unless 
the flux limiters are carefully constructed. When 
using a TVD scheme coupled with a multigrid 
acceleration it may be neccessary to limit the 
transfer of the residual to coarser meshes in the 
vicinity of shocks. Hence, there is a need for 
more work to extend the TVD theory to steady 


state problems and also weak shocks. 

• For time dependent flows we were not able to use 
a CFL greater than one even though the linear 
stability of the Runge-Kutta allowed larger time 
steps. For the steady state problems one can use 
larger time steps. Hence, it is efficient to use 
many stages in the Runge-Kutta method to in- 
crease the CFL stability limit even though one is 
not interested in high time accuracy. Neverthe- 
less, the limitations on the time step for time de- 
pendent problems indicates that even for steady 
state problems one should limit the local CFL 
near shocks to less than one at least in the tran- 
sient phase. This is crucial for hypersonic flow. 
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